home *** CD-ROM | disk | FTP | other *** search
/ Chip 2007 January, February, March & April / Chip-Cover-CD-2007-02.iso / Pakiet bezpieczenstwa / mini Pentoo LiveCD 2006.1 / mpentoo-2006.1.iso / livecd.squashfs / usr / lib / python2.4 / site-packages / Numeric / MLab.pyo (.txt) < prev    next >
Python Compiled Bytecode  |  2006-03-29  |  16KB  |  461 lines
  1. # Source Generated with Decompyle++
  2. # File: in.pyo (Python 2.4)
  3.  
  4. """Matlab(tm) compatibility functions.
  5.  
  6. This will hopefully become a complete set of the basic functions available in
  7. matlab.  The syntax is kept as close to the matlab syntax as possible.  One 
  8. fundamental change is that the first index in matlab varies the fastest (as in 
  9. FORTRAN).  That means that it will usually perform reductions over columns, 
  10. whereas with this object the most natural reductions are over rows.  It's perfectly
  11. possible to make this work the way it does in matlab if that's desired.
  12. """
  13. from Numeric import *
  14. import RandomArray
  15.  
  16. def rand(*args):
  17.     '''rand(d1,...,dn) returns a matrix of the given dimensions
  18.     which is initialized to random numbers from a uniform distribution
  19.     in the range [0,1).
  20.     '''
  21.     return RandomArray.random(args)
  22.  
  23.  
  24. def randn(*args):
  25.     '''u = randn(d0,d1,...,dn) returns zero-mean, unit-variance Gaussian
  26.     random numbers in an array of size (d0,d1,...,dn).'''
  27.     x1 = RandomArray.random(args)
  28.     x2 = RandomArray.random(args)
  29.     return sqrt(-2 * log(x1)) * cos(2 * pi * x2)
  30.  
  31.  
  32. def eye(N, M = None, k = 0, typecode = None):
  33.     '''eye(N, M=N, k=0, typecode=None) returns a N-by-M matrix where the 
  34.     k-th diagonal is all ones, and everything else is zeros.
  35.     '''
  36.     if M is None:
  37.         M = N
  38.     
  39.     if type(M) == type('d'):
  40.         typecode = M
  41.         M = N
  42.     
  43.     m = equal(subtract.outer(arange(N), arange(M)), -k)
  44.     return asarray(m, typecode = typecode)
  45.  
  46.  
  47. def tri(N, M = None, k = 0, typecode = None):
  48.     '''tri(N, M=N, k=0, typecode=None) returns a N-by-M matrix where all
  49.     the diagonals starting from lower left corner up to the k-th are all ones.
  50.     '''
  51.     if M is None:
  52.         M = N
  53.     
  54.     if type(M) == type('d'):
  55.         typecode = M
  56.         M = N
  57.     
  58.     m = greater_equal(subtract.outer(arange(N), arange(M)), -k)
  59.     return m.astype(typecode)
  60.  
  61.  
  62. def diag(v, k = 0):
  63.     '''diag(v,k=0) returns the k-th diagonal if v is a matrix or
  64.     returns a matrix with v as the k-th diagonal if v is a vector.
  65.     '''
  66.     v = asarray(v)
  67.     s = v.shape
  68.     if len(s) == 1:
  69.         n = s[0] + abs(k)
  70.         if k > 0:
  71.             v = concatenate((zeros(k, v.typecode()), v))
  72.         elif k < 0:
  73.             v = concatenate((v, zeros(-k, v.typecode())))
  74.         
  75.         return eye(n, k = k) * v
  76.     elif len(s) == 2:
  77.         v = add.reduce(eye(s[0], s[1], k = k) * v)
  78.         if k > 0:
  79.             return v[k:]
  80.         elif k < 0:
  81.             return v[:k]
  82.         else:
  83.             return v
  84.     else:
  85.         raise ValueError, 'Input must be 1- or 2-D.'
  86.  
  87.  
  88. def fliplr(m):
  89.     '''fliplr(m) returns a 2-D matrix m with the rows preserved and
  90.     columns flipped in the left/right direction.  Only works with 2-D
  91.     arrays.
  92.     '''
  93.     m = asarray(m)
  94.     if len(m.shape) != 2:
  95.         raise ValueError, 'Input must be 2-D.'
  96.     
  97.     return m[(:, ::-1)]
  98.  
  99.  
  100. def flipud(m):
  101.     '''flipud(m) returns a 2-D matrix with the columns preserved and
  102.     rows flipped in the up/down direction.  Only works with 2-D arrays.
  103.     '''
  104.     m = asarray(m)
  105.     if len(m.shape) != 2:
  106.         raise ValueError, 'Input must be 2-D.'
  107.     
  108.     return m[::-1]
  109.  
  110.  
  111. def rot90(m, k = 1):
  112.     '''rot90(m,k=1) returns the matrix found by rotating m by k*90 degrees
  113.     in the counterclockwise direction.
  114.     '''
  115.     m = asarray(m)
  116.     if len(m.shape) != 2:
  117.         raise ValueError, 'Input must be 2-D.'
  118.     
  119.     k = k % 4
  120.     if k == 0:
  121.         return m
  122.     elif k == 1:
  123.         return transpose(fliplr(m))
  124.     elif k == 2:
  125.         return fliplr(flipud(m))
  126.     elif k == 3:
  127.         return fliplr(transpose(m))
  128.     
  129.  
  130.  
  131. def tril(m, k = 0):
  132.     '''tril(m,k=0) returns the elements on and below the k-th diagonal of
  133.     m.  k=0 is the main diagonal, k > 0 is above and k < 0 is below the main
  134.     diagonal.
  135.     '''
  136.     svsp = m.spacesaver()
  137.     m = asarray(m, savespace = 1)
  138.     out = tri(m.shape[0], m.shape[1], k = k, typecode = m.typecode()) * m
  139.     out.savespace(svsp)
  140.     return out
  141.  
  142.  
  143. def triu(m, k = 0):
  144.     '''triu(m,k=0) returns the elements on and above the k-th diagonal of
  145.     m.  k=0 is the main diagonal, k > 0 is above and k < 0 is below the main
  146.     diagonal.
  147.     '''
  148.     svsp = m.spacesaver()
  149.     m = asarray(m, savespace = 1)
  150.     out = (1 - tri(m.shape[0], m.shape[1], k - 1, m.typecode())) * m
  151.     out.savespace(svsp)
  152.     return out
  153.  
  154.  
  155. def max(m, axis = 0):
  156.     '''max(m,axis=0) returns the maximum of m along dimension axis.
  157.     '''
  158.     m = asarray(m)
  159.     return maximum.reduce(m, axis)
  160.  
  161.  
  162. def min(m, axis = 0):
  163.     '''min(m,axis=0) returns the minimum of m along dimension axis.
  164.     '''
  165.     m = asarray(m)
  166.     return minimum.reduce(m, axis)
  167.  
  168.  
  169. def ptp(m, axis = 0):
  170.     '''ptp(m,axis=0) returns the maximum - minimum along the the given dimension
  171.     '''
  172.     m = asarray(m)
  173.     return max(m, axis) - min(m, axis)
  174.  
  175.  
  176. def mean(m, axis = 0):
  177.     '''mean(m,axis=0) returns the mean of m along the given dimension.
  178.        If m is of integer type, returns a floating point answer.
  179.     '''
  180.     m = asarray(m)
  181.     return add.reduce(m, axis) / float(m.shape[axis])
  182.  
  183.  
  184. def msort(m):
  185.     '''msort(m) returns a sort along the first dimension of m as in MATLAB.
  186.     '''
  187.     m = asarray(m)
  188.     return transpose(sort(transpose(m)))
  189.  
  190.  
  191. def median(m):
  192.     '''median(m) returns a median of m along the first dimension of m.
  193.     '''
  194.     sorted = msort(m)
  195.     if sorted.shape[0] % 2 == 1:
  196.         return sorted[int(sorted.shape[0] / 2)]
  197.     else:
  198.         sorted = msort(m)
  199.         index = sorted.shape[0] / 2
  200.         return (sorted[index - 1] + sorted[index]) / 2.0
  201.  
  202.  
  203. def std(m, axis = 0):
  204.     '''std(m,axis=0) returns the standard deviation along the given 
  205.     dimension of m.  The result is unbiased with division by N-1.
  206.     If m is of integer type returns a floating point answer.
  207.     '''
  208.     x = asarray(m)
  209.     n = float(x.shape[axis])
  210.     mx = asarray(mean(x, axis))
  211.     if axis < 0:
  212.         axis = len(x.shape) + axis
  213.     
  214.     mx.shape = mx.shape[:axis] + (1,) + mx.shape[axis:]
  215.     x = x - mx
  216.     return sqrt(add.reduce(x * x, axis) / (n - 1.0))
  217.  
  218.  
  219. def cumsum(m, axis = 0):
  220.     '''cumsum(m,axis=0) returns the cumulative sum of the elements along the
  221.     given dimension of m.
  222.     '''
  223.     m = asarray(m)
  224.     return add.accumulate(m, axis)
  225.  
  226.  
  227. def prod(m, axis = 0):
  228.     '''prod(m,axis=0) returns the product of the elements along the given
  229.     dimension of m.
  230.     '''
  231.     m = asarray(m)
  232.     return multiply.reduce(m, axis)
  233.  
  234.  
  235. def cumprod(m, axis = 0):
  236.     '''cumprod(m) returns the cumulative product of the elments along the
  237.     given dimension of m.
  238.     '''
  239.     m = asarray(m)
  240.     return multiply.accumulate(m, axis)
  241.  
  242.  
  243. def trapz(y, x = None, axis = -1):
  244.     '''trapz(y,x=None,axis=-1) integrates y along the given dimension of
  245.     the data array using the trapezoidal rule.
  246.     '''
  247.     y = asarray(y)
  248.     if x is None:
  249.         d = 1.0
  250.     else:
  251.         d = diff(x, axis = axis)
  252.     y = asarray(y)
  253.     nd = len(y.shape)
  254.     slice1 = [
  255.         slice(None)] * nd
  256.     slice2 = [
  257.         slice(None)] * nd
  258.     slice1[axis] = slice(1, None)
  259.     slice2[axis] = slice(None, -1)
  260.     return add.reduce(d * (y[slice1] + y[slice2]) / 2.0, axis)
  261.  
  262.  
  263. def diff(x, n = 1, axis = -1):
  264.     """diff(x,n=1,axis=-1) calculates the n'th difference along the axis specified.
  265.        Note that the result is one shorter in the axis'th dimension.
  266.        Returns x if n == 0. Raises ValueError if n < 0.
  267.     """
  268.     x = asarray(x)
  269.     nd = len(x.shape)
  270.     if nd == 0:
  271.         nd = 1
  272.     
  273.     if n < 0:
  274.         raise ValueError, 'MLab.diff, order argument negative.'
  275.     elif n == 0:
  276.         return x
  277.     elif n == 1:
  278.         slice1 = [
  279.             slice(None)] * nd
  280.         slice2 = [
  281.             slice(None)] * nd
  282.         slice1[axis] = slice(1, None)
  283.         slice2[axis] = slice(None, -1)
  284.         return x[slice1] - x[slice2]
  285.     else:
  286.         return diff(diff(x, 1, axis), n - 1)
  287.  
  288.  
  289. def cov(m, y = None, rowvar = 0, bias = 0):
  290.     '''Estimate the covariance matrix.
  291.  
  292.     If m is a vector, return the variance.  For matrices where each row
  293.     is an observation, and each column a variable, return the covariance
  294.     matrix.  Note that in this case diag(cov(m)) is a vector of
  295.     variances for each column.
  296.  
  297.     cov(m) is the same as cov(m, m)
  298.  
  299.     Normalization is by (N-1) where N is the number of observations
  300.     (unbiased estimate).  If bias is 1 then normalization is by N.
  301.  
  302.     If rowvar is zero, then each row is a variable with
  303.     observations in the columns.
  304.     '''
  305.     if y is None:
  306.         y = m
  307.     else:
  308.         y = y
  309.     if rowvar:
  310.         m = transpose(m)
  311.         y = transpose(y)
  312.     
  313.     if m.shape[0] == 1:
  314.         m = transpose(m)
  315.     
  316.     if y.shape[0] == 1:
  317.         y = transpose(y)
  318.     
  319.     N = m.shape[0]
  320.     if y.shape[0] != N:
  321.         raise ValueError, 'x and y must have the same number of observations.'
  322.     
  323.     m = m - mean(m, axis = 0)
  324.     y = y - mean(y, axis = 0)
  325.     if bias:
  326.         fact = N * 1.0
  327.     else:
  328.         fact = N - 1.0
  329.     val = squeeze(dot(transpose(m), conjugate(y)) / fact)
  330.     return val
  331.  
  332.  
  333. def corrcoef(x, y = None):
  334.     '''The correlation coefficients
  335.     '''
  336.     c = cov(x, y)
  337.     d = diag(c)
  338.     return c / sqrt(multiply.outer(d, d))
  339.  
  340. import LinearAlgebra
  341.  
  342. def squeeze(a):
  343.     '''squeeze(a) returns a with any ones from the shape of a removed'''
  344.     a = asarray(a)
  345.     b = asarray(a.shape)
  346.     return reshape(a, tuple(compress(not_equal(b, 1), b)))
  347.  
  348.  
  349. def kaiser(M, beta):
  350.     '''kaiser(M, beta) returns a Kaiser window of length M with shape parameter
  351.     beta. It depends on the cephes module for the modified bessel function i0.
  352.     '''
  353.     import cephes as cephes
  354.     n = arange(0, M)
  355.     alpha = (M - 1) / 2.0
  356.     return cephes.i0(beta * sqrt(1 - ((n - alpha) / alpha) ** 2.0)) / cephes.i0(beta)
  357.  
  358.  
  359. def blackman(M):
  360.     '''blackman(M) returns the M-point Blackman window.
  361.     '''
  362.     n = arange(0, M)
  363.     return (0.41999999999999998 - 0.5 * cos(2.0 * pi * n / (M - 1))) + 0.080000000000000002 * cos(4.0 * pi * n / (M - 1))
  364.  
  365.  
  366. def bartlett(M):
  367.     '''bartlett(M) returns the M-point Bartlett window.
  368.     '''
  369.     n = arange(0, M)
  370.     return where(less_equal(n, (M - 1) / 2.0), 2.0 * n / (M - 1), 2.0 - 2.0 * n / (M - 1))
  371.  
  372.  
  373. def hanning(M):
  374.     '''hanning(M) returns the M-point Hanning window.
  375.     '''
  376.     n = arange(0, M)
  377.     return 0.5 - 0.5 * cos(2.0 * pi * n / (M - 1))
  378.  
  379.  
  380. def hamming(M):
  381.     '''hamming(M) returns the M-point Hamming window.
  382.     '''
  383.     n = arange(0, M)
  384.     return 0.54000000000000004 - 0.46000000000000002 * cos(2.0 * pi * n / (M - 1))
  385.  
  386.  
  387. def sinc(x):
  388.     '''sinc(x) returns sin(pi*x)/(pi*x) at all points of array x.
  389.     '''
  390.     y = pi * where(x == 0, 9.9999999999999995e-21, x)
  391.     return sin(y) / y
  392.  
  393.  
  394. def eig(v):
  395.     '''[x,v] = eig(m) returns the eigenvalues of m in x and the corresponding
  396.     eigenvectors in the rows of v.
  397.     '''
  398.     return LinearAlgebra.eigenvectors(v)
  399.  
  400.  
  401. def svd(v):
  402.     '''[u,x,v] = svd(m) return the singular value decomposition of m.
  403.     '''
  404.     return LinearAlgebra.singular_value_decomposition(v)
  405.  
  406.  
  407. def angle(z):
  408.     '''phi = angle(z) return the angle of complex argument z.'''
  409.     z = asarray(z)
  410.     if z.typecode() in [
  411.         'D',
  412.         'F']:
  413.         zimag = z.imag
  414.         zreal = z.real
  415.     else:
  416.         zimag = 0
  417.         zreal = z
  418.     return arctan2(zimag, zreal)
  419.  
  420.  
  421. def roots(p):
  422.     ''' return the roots of the polynomial coefficients in p.
  423.  
  424.     The values in the rank-1 array p are coefficients of a polynomial.
  425.     If the length of p is n+1 then the polynomial is
  426.     p[0] * x**n + p[1] * x**(n-1) + ... + p[n-1]*x + p[n]
  427.     '''
  428.     if type(p) in [
  429.         types.IntType,
  430.         types.FloatType,
  431.         types.ComplexType]:
  432.         p = asarray([
  433.             p])
  434.     else:
  435.         p = asarray(p)
  436.     n = len(p)
  437.     if len(p.shape) != 1:
  438.         raise ValueError, 'Input must be a rank-1 array.'
  439.     
  440.     ind = 0
  441.     while p[ind] == 0:
  442.         ind = ind + 1
  443.     p = asarray(p[ind:])
  444.     N = len(p)
  445.     root = zeros((N - 1,), 'D')
  446.     ind = len(p)
  447.     while p[ind - 1] == 0:
  448.         ind = ind - 1
  449.     p = asarray(p[:ind])
  450.     N = len(p)
  451.     if N > 1:
  452.         A = diag(ones((N - 2,), p.typecode()), -1)
  453.         A[(0, :)] = -p[1:] / p[0]
  454.         root[:N - 1] = eig(A)[0]
  455.     
  456.     if (root.typecode() == 'F' or allclose(root.imag, 0, rtol = 9.9999999999999995e-08) or root.typecode() == 'D') and allclose(root.imag, 0, rtol = 1e-14):
  457.         root = root.real
  458.     
  459.     return root
  460.  
  461.